%{
For FBHALE software
Copyright (C) 2018-present Facebook, Inc.

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%}
  
function [ optim ] = GetSolarMassDistribution(optim)
%This function runs DiscretizeSolarPanels then computes the CG of chordwise
%strips of solar panels. It also computes the maximum input power on winter
%solstice to size the mppt. MPPT mass is distributed over the wing
%proportional to strip solar area the mppt mass is positioned to minimize the
%pitching moment generated by the solar panels. 

% Get Solar panel placement
optim = DiscretizeSolarPanels( optim );

fn = optim.solar.totalSolarfn ;
area = optim.solar.totalSolarAreaVect_m2;
Xb  = optim.solar.totalSolarXb ;

% Find furtherst forward and aft panel at each break point
% x = 0 is 1/4 chord for Xb 
rootX_m         =  Xb(Xb(:,2)  == 0,1);
breakPointX_m   =  Xb(Xb(:,2)  == optim.wing.bref_m/2*optim.wing.yobo2(2),1);
tipX_m          =  Xb(Xb(:,2)  == max(abs(Xb(:,2))),1);
solarLE = [max(rootX_m) max(breakPointX_m), max(tipX_m)];
solarTE = [min(rootX_m) min(breakPointX_m), min(tipX_m)];

% Shift to global axis sytem to LE = 0
solarLE = solarLE*-1;
solarTE = solarTE*-1;

% Break wing into 2 spanwise sections
sectionWidth_m  = optim.wing.bref_m/2/optim.wing.solar.nRepresentativePointMasses * ones(1,optim.wing.solar.nRepresentativePointMasses); 
pointMassY_m    = cumsum(sectionWidth_m)-sectionWidth_m(1)/2;

% Find section LE and TE edges
sectionSolarLE = interp1(optim.wing.bref_m/2*optim.wing.yobo2, solarLE, pointMassY_m);
sectionSolarTE = interp1(optim.wing.bref_m/2*optim.wing.yobo2, solarTE, pointMassY_m);

% Find panel CM
panelXCM_m    = (sectionSolarLE + sectionSolarTE)/2;
panelArea_m2   = ((abs(sectionSolarLE)+abs(sectionSolarTE)).*sectionWidth_m);
panelArea_m2   = panelArea_m2*(sum(optim.wing.solar.area_m2))/2/sum(panelArea_m2); %half span
panelMass     =  panelArea_m2*optim.solar.panelDensity_kg_m2;

%find mppt weight based on peak solar power on winter solstice
[optim.solar.mppt.mass_kg, optim.solar.mppt.mpptSizingPower_W] = GetMPPTMass(optim, fn, area);

% Find sectional mppt weight
mpptWeightPerSection = optim.solar.mppt.mass_kg/2/optim.wing.solar.nRepresentativePointMasses; % over 2 due to two wings

% Now attempt to cancel out moment created by the solar cells around the
% Elastic-axis using mppt placement
eaxis_x_m = -1*interp1(optim.wing.structure.sRefY_m{1} , zeros(1,length(optim.wing.structure.sRefX_m{1}))+optim.wing.structure.Cea_m{1}', pointMassY_m);

% Find moment generated by the solar panel about the s axis
solarPanelMoment = -1*(panelXCM_m-eaxis_x_m).*panelMass;

% Place mppt to counteract panel moment
mpptX_msaxisframe = -1*solarPanelMoment/mpptWeightPerSection;
furthestForwardMpptx_m = interp1(optim.wing.structure.sRefY_m{1}, (optim.wing.structure.Xax{1}-.05)'.*optim.wing.structure.chord_m{1}', ...
    pointMassY_m);
mpptX_msaxisframe(mpptX_msaxisframe>interp1(optim.wing.structure.sRefY_m{1}, ...
    (optim.wing.structure.Xax{1}-.05).*optim.wing.structure.chord_m{1}', pointMassY_m)) =  ... 
    furthestForwardMpptx_m(mpptX_msaxisframe>interp1(optim.wing.structure.sRefY_m{1}, ...
    (optim.wing.structure.Xax{1}-.05).*optim.wing.structure.chord_m{1}', pointMassY_m)) ;

% Find total mass and place to generate an equivalent moment
totalMoment = solarPanelMoment+mpptX_msaxisframe*mpptWeightPerSection;
totalSectionalMass_kg = panelMass+mpptWeightPerSection;

entities = optim.solar.surfacesWithSolar;
for e = 1:length(entities)
    for n = 1:length(optim.(entities{e}).N)
        % All values relative to 1/4 chord, actual postion set in
        % UpdateyLocations
        if strcmp(entities{e}, 'wing')
            optim.wing.solar.xRef_m     = -1*(totalMoment./totalSectionalMass_kg) + -1*eaxis_x_m; %relative to 1/4 chord
            optim.wing.solar.xRef_m   	= optim.wing.solar.xRef_m;
            optim.wing.solar.y_m     	= pointMassY_m;
            optim.wing.solar.mass_kg 	= totalSectionalMass_kg;
        else % htail and vtail
            optim.(entities{e}).solar.mass_kg   = sum(optim.(entities{e}).solar.area_m2)*optim.solar.panelDensity_kg_m2;
            optim.(entities{e}).solar.x_m       = 0; % relative to airfoil tail location
            optim.(entities{e}).solar.y_m       = 0;
        end

    end
end
end